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I. INTRODUCTION 

Structure formation in the Universe is a rich and fascinating problem, that touches many 
sides of Physics, from theories of the origin of primordial fluctuations, to detailed calcula- 
tions of galaxy dynamics in the present epoch, including radiation pressure from stars and 
absorption in clouds, star birth and star mass loss. It is of great topical interest today be- 
cause the recent and still improving observational data on the inhomogeneities in the cosmic 
microwave background radiation. 
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Structure formation is also a challenging problem of classical Physics, i.e. how a small 
perturbation of a spatially almost uniform universe develops, first linearly and then non- 
linearly, to the pronounced structures we see today. The theory of the linear stage of this 
process was developed by Lifshitz in the late 1940'ies, see |T8| , p!9| , while the nonlinear 
problem can either be treated as compressible hydrodynamics, or on the kinetic level by 
N-body calculations or Vlasov equations. In both always in non-linear hydrody- 

namics and gas dynamics of large systems with many scales, approximations have to be 
made. These typically amount to introducing turbulent viscosities or diffusivities, either 
directly in modeling, or indirectly in the numerical scheme. Even so, large simulations of 
nonlinear gravitational instabilities are numerically difficult problems, and a good deal of 
expertise and familiarity with a scheme is required to evaluate the status and validity of a 
result. 



For the reasons sketched above, there has since a long time been an interest in simpli- 
fied models of structure formation. The very simplest starts from the observation that in a 
one-dimensional setting, gravitational attraction is a Lagrangian invariant between particle 
collisions. This means that starting from an initially single-stream solution, where velocity 
is a function of position, one can straight-forwardly compute both the linear and the non- 
linear behavior up the moment of caustic formation, and deduce e.g. the spectra of density 
perturbations up to that time. After caustic formation several approaches are possible: ei- 
ther one may ignore the resulting change in the gravitational force (pancake, or Zeldovich 
model |jl6|), or one may assume that its effects may be modeled by an effective diffusive term 
(adhesion, or Gurbatov-Saichev-Shandarin model [^]). Both these models can be solved in 
one, two or three dimensions, and e.g. relations between statistics of initial conditions and 
statistics of the solutions computed. The last model, which gives Burgers' equation, leads 
already to interesting probabilistic and numerical problems, surveyed in |]T^ . 
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The present paper is the third in a series where we stay a httle closer to the original problem 
than working with Burgers' equation, and study directly the collision-less self-gravitating 
media. One motivation is to investigate if observed self-similar behavior of the solutions of 
Burgers' equation also appears in a self-gravitating system. In we introduced a simplified 
model for investigating the structure formation in an expanding Universe. This discussion 
will be taken up here again, in connection with the model of Rouet et al P, [l2| , p!3| for the 
growth of perturbations in a flat Universe. In particular, we will derive a new model that 



is also suitable for numerical application {Quintic model). In [0, in collaboration with 
A. Noullez, we presented an improved numerical algorithm, which takes 0{logN) opera- 
tions per collision to simulate particles. This new tool, extensively adopted here, allows 
a significant enhancement of statistics over our previous investigations. 



The paper is organized as follows. In section |T| we derive the Newtonian approximation, 
starting from General Relativity. In section |T| we consider the special case of a stratified 
perturbation in a 3D homogeneous expanding Universe. Then, in section ^V\, we re-derive 
the approximate solution presented in |]T|. In section |V| we specialize to a fiat Universe and 
discuss the model of P, p!^P^ , and the new one presented here. In section |V| we investigate 
properties of the mass distribution. Finally, in section [VII] we sum up and discuss our results. 

II. FROM GENERAL RELATIVITY TO CLASSICAL MECHANICS 

In General Relativity, points in space and time manifold are labeled by four coordinates 
x^. The Einstein conventions of covariant and contravariant indices and summation over re- 
peated indices are assumed. The field g^i, is (up to reparametrization) given by the solutions 
of Einstein's field equations 

Rf,u - ^9fiuR = ^-nGT^y (1) 

where Rfj,u is the curvature tensor, R = g^^^R^j^^ its trace, T^^ the energy-momentum tensor 
of the matter fields and G Newton's gravitational constant. 
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For a perfect fluid moving with four-velocity U relative to a given coordinate system, the 
energy momentum tensor is 

T>"' = pg^"" + {p + p)U^'U'' (2) 

where p and p are the pressure and the energy density. Einsteins's equations simplify if the 
Universe is assumed homogeneous and isotropic, and reduce then to Freidmann's equation 
for a scale factor a 

■2,1 SttG 2 /o^ 
a + k = -^pa^ (3) 

where k is the sign of the curvature of three-dimensional slices (—1, or 1). A spatially 
isotropic and homogeneous metric g is given by a, e.g. in the Robertson- Walker metric with 
gtt equal to one, grr = — — kr"^), gee = —a^r"^ and g^^ = — a^r^sin^6', all off-diagonal 
elements zero. In a curved Friedmann Universe [k = 1,-1) a{t) is its actual radius, but in 
a fiat Friedmann Universe {k = 0) only the ratios of a's at different times have meaning. 
Equation (|^) should be supplemented by the law of conservation of energy (pa^)' = —Spa^ 
and an equation of state, p = p{p) . 



If the energy density is dominated by non-relativistic matter, pressure can be neglected 
and p(t) = p(^o)(^^|^)~^- If the Universe is fiat {k = 0) Friedmann's equation has the simple 
solution a{t) = a{to){-^)^, where the reference time tQ, where p(to) and a(to) are given, 
satisfies to = (67rG'p(to))^^- The present belief is that the Universe is actually fiat or close 
to fiat. For the slightly more complicated expressions that hold in an over-critical or an 
under-critical Universe, see e.g. |T8|. 



Linear perturbations of the Einstein equations (for g), around the Friedmann solutions 
(where g is given by a and k), can be classified as scalar, vector and tensor, where the 
scalar perturbations are coupled modes of density and potential proper velocity, the vector 



modes are solenoidal proper velocity fields, and the tensor perturbations are gravitational 
waves. If we neglect gravitational waves, Friedmann's equations and the equations for the 



perturbations can in fact be derived in a Newtonian setting, as explained in |T8|, end of 
section 15.1, which we will briefiy recall here. 

Take a spherical part S of the Universe, which is assumed homogeneous and isotropic outside 
S. The solutions of Einstein's equations in S do then not depend on the value of the density 
outside S, just as there is no gravitational force acting in Newtonian gravity on the inside 
of a spherical shell. For small masses Einstein gravity is well approximated by Newtonian 
gravity. We first assume that mass density is constant in S, equal to Pb(t). If then S is 
sufficiently small and "r* is a coordinate system at rest with respect to the Universe outside 
S, and with origin in the center of S, a particle moves according to 

d'^T^ A-KmGr^ , , 1 , 



Introducing the proper coordinate it in S such that 

_^ a{t) 



a{to 

we can rewrite (HI) to 



(5) 



(fa AtiG , , 

^ = -^P.{t)a (6) 

Equation (P) is the time-time component of the Einstein field equations (|l|), assuming non- 
relativistic matter, and also the time derivative of (^, assuming mass conservation. 

We now turn to perturbations, and to simplify the notation we will write a for a(t), the 
scale factor at time t, and ao and po for a(to) and pbito), the scale factor and the background 
density at the reference time to- Let the density in S at time t fiuctuate around the mean 
value pb{t)- The Newtonian equations of motion for particles follow from a Lagrangian 

£ = ^ r?^ - mi(/)( r?,t) (7) 



where = inGp. By substituting the proper coordinate if and applying the canonical 
transformations 



we have 



X ^ 1 • 2 

£ = 2_^-m^Xi -mLp{xi) (9) 



where 99 = + |ff 3^^- The field equation for the new potential cp reads: 



Vl^ = AnGp^ + ^, (10) 



which, recalling (P), can be written |TT[| : 



a 
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V^¥. = 47rG - (p(ir,t)-p,(t)). (11) 
veto/ 

Equation ( p^T]) is the starting point of our work. The source of if is not the density itself, 
but the density perturbation {p~ pb{t)). This is as should be: in the absence of density 
fiuctuations each particle moves according to the uniform expansion of the Universe as a 
whole, but remains undisturbed in the proper coordinate if. 

The equation of motion associated to the Lagrangian (^ are 

i (^^] + ^M^) = 0, (12) 



dt \al / 

where, for simplicity, we neglect from here on the label i. The peculiar velocity of a particle 
is not if but if = -^if- The equations of motion in Eulerian form, for the peculiar velocity 
if, are therefore 



Equations (|T3]) (or ([T^)), (|TT]), and the equation of continuity define the Newtonian model 
of structure formation we study. Closing this section, it is appropriate to summarize in 
what respects the model is inexact: i) it ignores gravitational waves, which are degrees of 
freedom of the Einstein equations with no counter-parts in Newtonian theory; ii) it is limited 
to weak gravitational fields, in the sense that Newtonian mechanics holds for the peculiar 
velocity; iii) it does not treat correctly the expansion of a local patch with a density different 
from the average, since density perturbations act as sources of Newtonian (not Einsteinian) 
gravity. Points i) and ii) are not serious for our purposes. Point iii) is actually not as 
problematic as it seems, since the linear growth rates in the Newtonian model agrees with 
the growth rates in Lifshitz' full theory of linear perturbations to the Friedmann solutions of 
the Einstein equations, so would only be pertinent for the nonlinear growth of perturbations 
of sufficiently long wave length. 



We now consider the special case of a stratified perturbation: the velocity has one com- 
ponent only and varies with respect to this direction. Initial density also only varies in this 
direction. In the point particle picture, the density profile which goes into the Lagrangian 



where x is the comoving coordinate, in the direction of which the density and velocity varies. 
The Poisson equation leads to the following expression for the gravitational potential: 



III. STRATIFIED PERTURBATIONS 



(0) is 




(14) 




(15) 
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The integral in ([I5D should be taken over an interval from —L/2 to L eventually taken 
to infinity, or to where the perturbation vanishes Equation (|12|) reads 



dp X Qj dx (2 \ 

— + 2-— - 4'KGpb{t)x = - Egravix, t) (16) 

at'' a at \ao J 

where 

Egrav{x, t) = —2nG mjsign(x — Xj). (17) 
j 

The interesting thing is now that by the equation of continuity 

the time-dependences of AixGphit) and Egrav{x,t) are actually the same. As pointed 

out by Rouet et al ll^Jl^) time-dependence can then be concentrated in the term 2-"^^ 



' a dt 

by a suitable nonlinear transformation of the time variable. The choice should be 

* = (^) * 

where r has dimension of time. Equation (|I6|) is then transformed to 

d'^x dJU dx 

+ ~ 47rGpoa; = Eg^av (x, r) (20) 

The interest of this formulation is that, as in classical self-gravitating systems in one dimen- 
sion, Egyav is a Lagrangian invariant, proportional to the net mass difference to the right 
and to the left of a given particle at a given time. Therefore, as far the particle do not 
experience any crossing, the equation of motion (|20| ) reduces to the compact form: 



— + B(t)— - Cx = const, (21) 
ctr^ dr 

where B{t) = and C = AnGpQ. We have hence reduced the motion between collisions 
to a Ricatti equation. When that can be integrated in an efficient manner, we can therefore 
solve general Newtonian model with a fast event-driven scheme, as in ,|r^,|T^] . We will 
consider special cases when this can be done below. 
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IV. APPROXIMATE SOLUTION: FL MODEL 



In our previous paper we introduced a model we here call Friction Less (FL) model. 
That model basically cuts short the preceding discussion, and starts directly from equa- 



tion (12) with the assumption that we consider only phenomena on time scales much less 
than the age of the Universe. This means that we can take a constant, alternatively that 
our time scales must be smaller than a/ a. By a change of the spatial scale we can set a to 
one, and equation (|IB]) takes the much simpler form 

— - AirGpox = Egrav{x, t) FL model (22) 

Between collisions, (^) has solutions that are linear combination of two exponentials with 



rates ±A/47rGpo- This simple form allows the implementation of a fast and exact, up to 
round-off errors, integration scheme |jl|J^: the evolution of the system is mapped from one 
crossing to the next and the times the events occur are analytically computed by solving a 
quadratic equation. In the present paper we use a version of the new heap-based algorithm 
introduced in |T^. 



In our numerical experiment we consider a periodic perturbation with size 2L with re- 
flexion symmetry. Therefore we can simulate half of the symmetric system, confining N 
particles in a box with reflecting edges in —L/2 and L/2. The mass of the particles is 
m = N~'^. We choose as unit of length the spatial interval in which the particles are initially 
distributed, and, thus, the initial density po is equal to one. A natural choice of time scale 
is uj^ = (47rGpo)^^^^ = -y/i^o, where ujo is the Jeans frequency. 

In Fig. |l] we represent the evolution of the system in the phase space, starting with spatial 
uniform distribution. The initial velocity is a smooth function of position (first plots in Fig. 
|T]). As the time increases the multi-stream solution takes place. The spiral-like behavior 
develops and the thickness of the region where velocity is not single-valued grows quickly. 
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FIG. 1. Phase space portraits for the FL model with = 200. The plots refer to different t, 
expressed in unities of Jeans' time. The first plot shows the initial velocity profile. Positions and 
velocities are in arbitrary units. 



A. Zeldovich approximation and FL model 

As a test of the validity of the FL model we establish a connection with the Zeldovich 
approximation, closely parallel to the similar relation for self-gravitating systems without 
expansion and Burgers' equation. For an extensive recent review on the subject the reader 



can refer to [jl^. Consider the evolution of the system before the first crossing. According 



to eq. (P2D, we have: 



— = uj%{x~xo), (23) 



dt 

where Xq = x{0). The solution of eq. ( pSj) can be written in the form: 
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X = Xq + smh.{ujjt) (24) 

with vq = v{0). Following a suggestion by S.N. Gurbatov, we then introduce the new time 
defined as 6* = and the rescaled velocity U = v/6 = v/ cosh(u;ojit) [0. Thus, we get: 



X = Xq + U {xo)9 

(25) 

U{t,xo) = U{xo), 
which can be combined to give the Burgers equation: 

Thus, the evolution of a particle can be thought as a free motion, by performing an appro- 
priate rescaling of time and velocity, as far as no crossing takes place. As a less obvious 
consequence, the dynamics of the system can be described by macro-particles, carrying a 
significant part of the whole mass of the system. These macro-particles are initially located 
at the points of future caustic formation and evolve according to (^), as far they do not 
cross each other. Thus, we have access to information on the times of merging of structures 
without investigating the details of their inner dynamics. 

As an application consider the case of the smooth initial condition represented in Fig. |^. 
The diamonds states for two macro-particles of mass Nm/2. After a time = 3.3u!Jq the 
macro-particles have collided, and, as expected, the two multi-streamed structures have 
merged together. Time elapsed in original coordinate is t = 1.9ujJq. 
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FIG. 2. The first plot shows the initial condition in the phase space. The two diamond are the 
macro-particles of mass initially located in correspondence of the regions where the first 

singularities takes place. After 6 = 3.3wJq^, the macro-particles have collided, and, as predicted, 
the multi-stream structures merge {t = ujjQsinh^^{ujjQ6) = 1.9uJjq). Positions and velocities are 
in arbitrary units. 



V. MODELS OF A FLAT EXPANDING UNIVERSE 

In the matter-dominated era (pressure negligible) in a flat (critical) universe (fc = 0), the 
scale factor a(t) grows with time as a power-law 



^,2/3 



a{t) = ao (^-J . (27) 

Here to = (SvrGpo)^^^^, as follows from equation (H). Therefore, as also pointed out by Rouet 
et al |T2|JT3[] , the sole remaining time-dependent term in ( pOD becomes time- independent. Let 
us note that the substitution (p!9| ) can be integrated to 

r = tolog^ t = toexp(f) (28) 
12 



if we make the natural non-restrictive choice that T{to) = 0. Equation (pO]) then takes the 
form: 

(fx 1 dx 2 rn / \ ^ 1 1 / \ 



We refer to (p9|) as to the Quintic (Q) model. Between coUisions the right-hand side of (E9| 



is constant and the homogeneous equation has solutions 

x(r) = ci exp(^^^-^) + C2 exp(-^^^) (30) 

where ci and C2 are determined by x{tq) and x(to) in the transformed coordinate. The form 
of equation (PO) suggests introducing an auxiliary variable z = exp( ^^~^°-' ). The crossing 



3to 

times between consecutive particles are hence computed by solving numerically a quintic 
equation in the form: 

biz' - b2Z^ + h = 0, (31) 
where the coefficients bi, 62, ^3 are fixed by the states of the particles at the time of the last 



crossing. Therefore, the event-driven scheme |10] can be adopted to follow the dynamics of 



this system. In comparison, however, the FL model discussed above obviously leads to a 
quadratic equation, and the RF model presented below, which has been the benchmark in 
the literature, leads to a cubic equation. Both of these can be solved algebraically, while 
the quintic cannot. The detail of the numerical implementation are given elsewhere 0]. In 
the present numerical experiments we consider a system of N particles, all of the same mass 
m = N^'^, confined in a finite box with refiective boundaries conditions. The unit of length 
is the box size and time is expressed in unit of the inverse of the Jeans frequency ujjq. In 
Fig. ^ the evolution of the system in phase space is represented. The initial condition is the 
one we used for Fig. |l[ As is well-known, the time evolution leads to massive central core, 
and the system displays the expected spiral shape. Nevertheless, due to the action of the 
friction term, the particles remain well concentrated in the inner zone and the thickness of 
the multi-stream grows much slower than in the FL. 
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FIG. 3. Phase space portraits for Q model with N = 200. The plots refer to different r, 
specified in the figure legends in unities of Jeans' time. In terms of cosmological time, the times 
elapsed after to are respectively S.TOcjJq, 12.08ljJq and AQAJloJq. The first plot shows the initial 
velocity profile at to- Positions and velocities are in arbitrary units. 

Let us note the fact, that if we re-express in the original time variable t using 
and assume small initial perturbations of a uniform state, we then have, before any collisions, 

.(i).a.(i)''%a.(i)". (32) 

By the usual Eulerian-Lagrangian transformation, this means that small density perturba- 
tions develop as 

5p{x,t) = pi{x,to) (^^^ + p2{x,to) (33) 
Similarly, small peculiar velocities transform according to 

V{x,t) = Viix,to) i-j +V2{x,to)i-] (34) 
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In ( |33D and (Q) {pi,vi) is the growing mode of coupled density and velocity perturbations, 
while {p2,V2) is the decaying mode. It is easily checked that the growth rates in (^) and 
(p^) agree with the ones computed directly on the Eulerian side, see [|18[ or p9[] . 



A. The Rouet-Feix Model (RF) 

A very interesting model was introduced about ten years ago by Rouet, Feix and col- 
laborators |jl2|,|l3|, see also the recent paper 0]. One way to understand this model, in the 



spirit of the original derivation, is in the sense of a strictly Newtonian model of an expanding 
universe. Indeed, one may then specify an expansion 

/ t \^/^ 

«W = ao i^-^j , (35) 

without necessarily demanding, since one does not require compatibility with general rela- 
tivity, that the reference time is connected to the density like in to = {6nGp{to))~^^'^- 
We label the time t^^ since all change we actually have to do is to substitute (^ by the 
analogous ID equation 

^ = -AnGpit)a. (36) 

From (|35|) and (^6]) follow then 

p = {ISrrGfr' C = {ISnGpo)-'/' = ^u;J^ (37) 
where cujo = (AnGpoY^'^ Q Equation ([T9| ) can be integrated to 

r = t]^ log ^ = -^to log (^Y^ ' 



(38) 



^ Note that the numerical coefficient linking tj^ to Jeans' time uj j} is equal to the inverse of the 



numerical constant a introduced by the authors in the derivation in |13]. 
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which is hnearly related to the r given by (pHI). Between colhsions equation ( ]20|) is trans- 
formed into: 

d'^X UJjodx 2 . T3T7 J 1 /on\ 

-— r- H ^- uj jnX = const. Kb model, (39) 



which is the equation in [|T210. The solution of the homogeneous equation takes the form: 



x{t) = ciexp{^—-^—^) + C2exp(-A/2(r - to)ujjo), (40) 

which suggests the introduction of z = exp( '''^~^'^'^° ), in terms of which one should solve 
the cubic equation 

biz"^ - b2Z^ + &3 = (41) 

to find the particle crossings. As in the previous case, the coefficients hi, &2, &3 are fixed by 
the states of the particles at the time of the last crossing. We now transform back to time 
the coordinate t, using the second equahty of (pH), and find 



X 1/3 / N -2/3 

x{t) = c,{^] . (42) 



The exponents for the growing and decaying modes are different than in the Lifshitz-Bonner 
theory. Hence, the Rouet-Feix model does not predict correctly linear growth of small 
perturbations, and can therefore not be taken seriously as a quantitative model structure 
formation in a pressure-less gas in an expanding universe, obeying general relativity. 

Nevertheless, the RF model has the important merit of producing a final tractable ex- 
pression for computing algebraically the particle crossing times, and has qualitatively all the 
right ingredients of the dynamics (both friction and background terms). It can therefore be 



simulated by a straight-forward generalization of our heap-based event-driven scheme []T0 
For a more detailed discussion of the implementation, see [§]. 



We note now, that at least in linear theory, there is another sense one can give to 
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the RF model. Namely, if g is a wave-vector of a perturbation in a comoving frame 
{5p ~ Pgexp(?xg)), and we assume a non-negligible pressure p ~ p'^ with 7 = |, then 
such perturbations grow as 



p r^t a = --±{—- A ) ' A = — [4:6) 

o 60 



where Vg is the sound speed, equal to y-^- These exponents can be matched with RF if 
A = 2/3. Hence, the RF model describes the linear growth of perturbations of physical 



wave length ^tvg, which is \/S/2 times larger than the instantaneous Jeans length. We note 
furthermore that the Jeans' length, Ij = Vgj y/AnUp, scales with time as t"^^^ under the as- 
sumptions above, and that therefore the wave length of a perturbation remains in constant 
proportion to the Jeans wave length. The RF model may hence well have some validity 
also beyond linear theory. Let us add that the modifications of the values of the exponents, 
of (|4^ ) compared to (^2]), are quite reasonable. As the length scale of a perturbation ap- 
proaches the Jeans length from above, the difference between expanding and decaying rates 
diminishes, to eventually give rise to purely oscillatory behaviour below the Jeans' length. 



We now turn to simulations of the RF model. We consider a system of particles, all 
of the same mass m = A^"^, confined in a finite box with reflective boundaries conditions. 
The unit of length is the box size and time is expressed in unit of the inverse of the Jeans 
frequency ujq. In Fig. ^ we show the evolution of the system in the phase space, starting 
with the same initial condition we used for Fig. |l|. The result qualitatively agrees with the 
plot of Fig. ^ As for the Q model, due to the presence of the friction, the system develops 
a multi-stream region, that is more narrow than the one displayed by the FL. Nevertheless 
it has to be noticed that the process of structure formation is considerably slower, for the 
RF, than for the Q model. Similar phase space portraits are displayed at equal rescaled 
r, but within the framework of the Q model, shorter cosmological times t are required to 
attain such stages (see captions of Figs. |^, ^). In addition we observe that, for similar r, 
the thickness of the agglomeration is smaller in Fig ^ than in Fig. ^. A simple explanation 
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is, indeed, provided by looking to the relative importance of the friction and background 
terms, respectively in Q and RF models. In the latter, in fact, the role of the friction is 
overestimated and, as a consequence, the process of particles sticking is more pronounced. 
This observation suggests that, possibly, the adhesion model would result more similar to 
the approximate RF model, than to the exact Q model [Q. 
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FIG. 4. Phase space portraits for RF model with N = 200. The plots refer to different r, 
specified in the figure legends in unities of Jeans' time. In terms of cosmological time, the times 
elapsed after to are respectively 21.4632u;J'(|, 61. 9916a; Jq^ and 517.1AuJjq. The first plot shows the 
initial velocity profile at to. Positions and velocities are in arbitrary units. 



VI. MASS OCTAVE FUNCTION 

In this section we study the mass distributions for the FL, Q and RF models. Simulations 
are performed taking particles initially uniformly distributed on a line of size L. Thus, 
the initial inter-particles distance is Aa; = L/N. Velocities are generated as a Brownian 
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random process. This is done in the Fourier space representation where: 

v{x) = T.^Vke'"'^ (44) 

The sum over k extends from — vr / Aa; to vr / Ax in steps of 2 tt / L and v^k = v^. The Fourier 
components of positive k are then chosen as a random Gaussian independent variables with 
variances: 

where h is the Hurst exponent, the scahng exponent of the second order structure functions. 
The choice h = 0.5 corresponds to standard Brownian motion. The field generated by (pif ) 



and p5|) will be periodic with period L. We scale L to be one, and consider, as above. 



reflecting boundaries. 



Fig. 1^ shows a phase-space portrait for the FL model (t = A.OuJ^). Spiral structures 
are displayed, as well as stretched filaments connecting the dense agglomerations. Fig. |^ 
represents the phase-space for the Q model at r = A.OcUj^: the large-scale structures recall 
here even more the ones obtained in the framework of the adhesion model, even if it should 
be pointed out that, locally, multi-stream behaviour has in fact already started to develop 
(see inset in Fig. Substantially similar, from a merely qualitative point of view, is the 
phase space portrait for the RF model, at r = A.OcuJ^ (see Fig. 



To investigate quantitatively the particle distribution, we introduce the mass octave func- 
tion (MOF), already discussed in MOF measures the probability to find a non-zero 
contribution to the mass density, as function of the mass itself, coarse-grained in octaves. 
In other words, P(Am) is the cumulative probability of finding Am in an octave interval 
[mk,mk+i[, where rrik = {k = 0, 1,...). Practically, the MOF analysis reveals the degree 
of dishomogeneity of a probability distribution, since in terms of the MOF a uniform dis- 
tribution would correspond to a logarithmic histogram with only one non-zero entry. The 

19 



main plot of Fig. ^ show, in doubly logarithmic scale, the MOF computed for RF at an 
intermediate stage of evolution (r = S.SuJq). In a finite range the mass octave function 
displays a power-law decay with exponent —0.5. This result is in agreement with PJT^ 
where, in the framework of the adhesion model, the number density per unit length of shock 
locations holding mass m is shown to be distributed as power-law m~^~^. 



Figs |] show the MOF for a later time evolution of the RF model. An approximate power 
law regime is found (see small inset), and a numerical fit gives exponent a^^ = —0.25. 

The main plot of Fig |1^ represents the MOF for the Q model at rescaled time r = S.OujJq. 
There is no evidence of a transient power-law state, with slope —0.5, as shown for the RF. 
This observation again indicates that the adhesion picture is, in fact, closer to the RF model, 
than to the exact Q approach. The small inset displays, in doubly logarithmic scale, the 
MOF for the Q model at r = 3.0ujJq and r = 5.5uJq. In the latter case an approximate 
power-law, over a finite mass support, is displayed and the best numerical fit gives exponent 
= 0.1. 



Figs. |TT| shows on the other hand the MOF for the FL model. At an intermediate stage, 
(r = A.OluJq) the MOF is clearly not uniform, but best described as a structure with two 
peaks. At later times an approximate power-law regime is found, with a positive exponent 



{a^L ^ 0.4, at r = 5.5ujj^). 



We observe that the exponent is practically equidistant from the values of a and 
a^^. Nevertheless is positive, and therefore, concerning the mass aggregation, the exact 
Q dynamics is shown to be qualitatively more similar to the FL model, than to the RF. This 
is, indeed, a very surprising conclusion. 
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FIG. 5. Velocity field vs. position, for FL model, starting from a single-speed Brownian motion 
initial conditions (small inset). Here N = 16384 and t = A.OlOjq. Reflecting boundaries are 
assumed. Positions and velocities are in arbitrary units. 
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FIG. 6. Velocity field vs. position, for Q model, starting from a single-speed Brownian motion 
initial conditions. Here N = 16384 and r = 4:.0u}'Jq (or t = 109. 53a; Jq^). Reflecting boundaries are 
assumed. The inset is a zoom of a massive agglomeration. Positions and velocities are in arbitrary 
units. 
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FIG. 7. Velocity field vs. position, for RF model, starting from a single-speed Brownian motion 
initial conditions. Here N = 16384 and r = A.OcoJq (or t = 2.283 x IO^Wjq). Reflecting boundaries 
are assumed. The inset is a zoom of a massive agglomeration. Positions and velocities are in 
arbitrary units. 
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FIG. 8. Log-log normalized mass distribution in octaves (MOF) for Brownian initial velocities 
for RF (dashed line) at time r = S.SlOjq (or t = 517.14a; The number of particles is 16384, 
number of independent realizations 1000. The height of a column is the fraction of a total number 
of bins containing mass in the range [m,2m\. The bin size is I = 0.00195, a uniform density 
would hence correspond to a single column at abscissa 0.00195 and ordinate 1.0. An intermediate 
power-law regime is clearly displayed for RF model. The best numerical fit gives exponent a = —0.5 
(solid line). The small inset represents the MOFs in log-lin scale. Mass is in arbitrary units. 
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FIG. 9. Normalized mass distribution in octaves (MOF) for Brownian initial velocities at 
T = 5.5ujjQ (or t = 5.5007 x IO^wJq^), for RF model. The number of particles is 8192, number of 
independent realizations 1000. The bin size is Z = 0.01562. The small inset is a plot of the MOF 
in log-log scale. A numerical fit to a power law gives exponent a^^ = —0.25. Mass is in arbitrary 
units. 
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FIG. 10. Normalized mass distribution in octaves (MOF) for Brownian initial velocities at 
T = 3.0a; jQ (or t = 32.18a;jQ^), for Q model. The number of particles is 8192, number of independent 
realizations 500. The bin size is Z = 0.003906. The small inset represents the MOF in log-log scale 
for respectively t = 3.0a; Jq^ (dashed line and triangles) and r = S.bujjQ {t = 687.71a; Jq^, solid line 
and circles) . Bin size, number of particles and number of independent realizations as for the main 
figure. An intermediate power-law regime, with positive exponent, develops as the evolution is 
pushed forward. The best numerical fit gives exponent = 0.1 at r = b.bcUjQ . Mass is in 
arbitrary units. 
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FIG. 11. Normalized mass distribution in octaves (MOF) for Brownian initial velocities at 
t = 4:.0ujJq, for FL model. Number of particles and number of realizations as in Fig. ^ The bin 
size is / = 0.01562. The small inset represents the MOF in log-log scale for respectively t = A.OuJjq 
(dashed line and triangles) and t = S.ScjJq (solid line and circles). Bin size, number of particles 
and number of independent realizations as for the main figure. An intermediate power-law regime 
develops as the evolution is pushed forward. The best numerical fit gives exponent a^^ = 0.4 at 
t = 5.5ujj} . Mass is in arbitrary units. 



VII. CONCLUSIONS 

In this paper we have discussed the problem of structure formation in a three-dimensional 
expanding Universe, focusing on stratified perturbations in a collision-less medium. We de- 
rived a representation of the Lagrangian development of such perturbations from the linear 
regime to large times, which forms the basis of an efficient numerical scheme {Q model). In 
practice, this scheme requires solving automatically a large number of quintic equations, the 
coefficients of which are only known at run-time, and therefore entails algorithmic problems 
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that are discussed in [Q. 



We then re-derived as approximate solutions the Friction Less Model (FL), recently intro- 
duced in [|I|, and the Rouet-Feix (RF) model, and compared them to the exact Q dynamics. 
Quantitatively, the growth of small perturbations differs between RF and Q; the RF model 
is in this sense more similar to growth of finite wave length perturbations in a model with 
non-zero pressure. Numerical simulations have been carried out for smooth and random 
initial perturbations, and phase space portraits, respectively for FL, RF and Q, have been 
computed. The latters show compact collapsed structures, visually more like the shock 
waves (singular mass agglomerations) in the adhesion model (Burgers' equation). 



To make our analysis more quantitative, we studied the statistics of mass distributions from 
random initial conditions. For intermediate times, when the multi-stream behaviour is not 
well developed, the RF model agrees well with the scaling law for the frequency of agglom- 



erations of given size containing mass m in Burgers equation |TB|,|T^,|T^ , i.e. 



Sim) ~ m 



-l-h 



This is not true for either the Q or the FL model. In particular, the FL seems to tend to a 
bifractal structure, with either very high or very low density agglomerations. 



We finally studied the late time behaviour of the FL, Q and RF models. All then dis- 
play states which seem to have approximate power-law MOF. The actual states are however 
different, with the most numerous mass agglomerations in FL (positive a^^) having about 
equal mass (they then also carry most of the mass), while the most numerous mass agglom- 
erations in the RF (negative a^^) are still small, and carry a small proportion of the mass. 
The exponent a^, characteristic of the Q dynamics, is practically equidistant from a^^ and 
a^^. Nevertheless a*^ is positive and, therefore, surprisingly enough, we are led to conclude 
that, concerning the mass aggregation, the exact Q dynamics is qualitatively more similar 
to the FL model than to the RF model. 
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Summing up, we have compared three models for studying the evolution of a stratified 
perturbation in an expanding universe. Only the Q model is in quantitative agreement 
with the linear theory of growth of perturbations to the Friedmann solutions of the Einstein 
equations. Quantitatively, the approximate RF is shown to agree with the adhesion model, 
with respect to the mass aggregation. 
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